2D decaying turbulence
This example can be run online via . Also, it can be viewed as a Jupyter notebook via
.
A simulation of decaying two-dimensional turbulence.
using FourierFlows, Printf, Random, Plots
using Random: seed!
using FFTW: rfft, irfft
import GeophysicalFlows.TwoDNavierStokes
import GeophysicalFlows.TwoDNavierStokes: energy, enstrophy
import GeophysicalFlows: peakedisotropicspectrumChoosing a device: CPU or GPU
dev = CPU() # Device (CPU/GPU)Numerical, domain, and simulation parameters
First, we pick some numerical and physical parameters for our model.
n, L = 128, 2π # grid resolution and domain length
# Then we pick the time-stepper parameters
dt = 1e-2 # timestep
nsteps = 4000 # total number of steps
nsubs = 20 # number of steps between each plotProblem setup
We initialize a Problem by providing a set of keyword arguments. The stepper keyword defines the time-stepper to be used.
prob = TwoDNavierStokes.Problem(dev; nx=n, Lx=L, ny=n, Ly=L, dt=dt, stepper="FilteredRK4")Next we define some shortcuts for convenience.
sol, clock, vars, grid = prob.sol, prob.clock, prob.vars, prob.grid
x, y = grid.x, grid.ySetting initial conditions
Our initial condition closely tries to reproduce the initial condition used in the paper by McWilliams (JFM, 1984)
seed!(1234)
k₀, E₀ = 6, 0.5
ζ₀ = peakedisotropicspectrum(grid, k₀, E₀, mask=prob.timestepper.filter)
TwoDNavierStokes.set_zeta!(prob, ζ₀)Let's plot the initial vorticity field:
heatmap(x, y, vars.zeta',
aspectratio = 1,
c = :balance,
clim = (-40, 40),
xlims = (-L/2, L/2),
ylims = (-L/2, L/2),
xticks = -3:3,
yticks = -3:3,
xlabel = "x",
ylabel = "y",
title = "initial vorticity",
framestyle = :box)
Diagnostics
Create Diagnostics – energy and enstrophy functions are imported at the top.
E = Diagnostic(energy, prob; nsteps=nsteps)
Z = Diagnostic(enstrophy, prob; nsteps=nsteps)
diags = [E, Z] # A list of Diagnostics types passed to "stepforward!" will be updated every timestep.Output
We choose folder for outputing .jld2 files and snapshots (.png files).
filepath = "."
plotpath = "./plots_decayingTwoDNavierStokes"
plotname = "snapshots"
filename = joinpath(filepath, "decayingTwoDNavierStokes.jld2")Do some basic file management
if isfile(filename); rm(filename); end
if !isdir(plotpath); mkdir(plotpath); endAnd then create Output
get_sol(prob) = prob.sol # extracts the Fourier-transformed solution
get_u(prob) = irfft(im * grid.l .* grid.invKrsq .* sol, grid.nx)
out = Output(prob, filename, (:sol, get_sol), (:u, get_u))
saveproblem(out)Visualizing the simulation
We initialize a plot with the vorticity field and the time-series of energy and enstrophy diagnostics.
p1 = heatmap(x, y, vars.zeta',
aspectratio = 1,
c = :balance,
clim = (-40, 40),
xlims = (-L/2, L/2),
ylims = (-L/2, L/2),
xticks = -3:3,
yticks = -3:3,
xlabel = "x",
ylabel = "y",
title = "vorticity, t=" * @sprintf("%.2f", clock.t),
framestyle = :box)
p2 = plot(2, # this means "a plot with two series"
label = ["energy E(t)/E(0)" "enstrophy Z(t)/Z(0)"],
legend = :right,
linewidth = 2,
alpha = 0.7,
xlabel = "t",
xlims = (0, 41),
ylims = (0, 1.1))
l = @layout Plots.grid(1, 2)
p = plot(p1, p2, layout = l, size = (900, 400))
Time-stepping the Problem forward
We time-step the Problem forward in time.
startwalltime = time()
anim = @animate for j = 0:Int(nsteps/nsubs)
if j % (1000 / nsubs) == 0
cfl = clock.dt * maximum([maximum(vars.u) / grid.dx, maximum(vars.v) / grid.dy])
log = @sprintf("step: %04d, t: %d, cfl: %.2f, ΔE: %.4f, ΔZ: %.4f, walltime: %.2f min",
clock.step, clock.t, cfl, E.data[E.i]/E.data[1], Z.data[Z.i]/Z.data[1], (time()-startwalltime)/60)
println(log)
end
p[1][1][:z] = vars.zeta
p[1][:title] = "vorticity, t=" * @sprintf("%.2f", clock.t)
push!(p[2][1], E.t[E.i], E.data[E.i]/E.data[1])
push!(p[2][2], Z.t[Z.i], Z.data[Z.i]/Z.data[1])
stepforward!(prob, diags, nsubs)
TwoDNavierStokes.updatevars!(prob)
end
gif(anim, "twodturb.gif", fps=18)